In-silico selection of peptides for the recognition of imidacloprid

The sensitive detection of pesticides using low-cost receptors designed from peptides can widen their uses in the environmental surveillance for emerging pollutants. In-silico selection of peptides can help accelerate the design of receptor sequence banks for a given target of interest. In this work, we started from Lymnaea stagnalis acetylcholine-binding protein Q55R mutant receptor-imidacloprid complex, available in the PDB databank, to select three primary short peptides (YSP09, DMR12, WQW13 respectively having 9, 12 and 13 amino acids (AA) in length) from the pesticide interacting zones with the A, B and C chains of the nicotinic receptor. Using molecular docking and molecular dynamics (MD) simulations, we showed that the three peptides can form complexes with the target imidacloprid, having energies close to that obtained from a reference RNR12 peptide. Combination of these peptides allowed preparing a new set of longer peptides (YSM21, PSM22, PSW31 and WQA34) that have higher stability and affinity as shown by the MM-PBSA calculations. In particular, the WQA34 peptide displayed an average binding free energy of –6.44±0.27 kcal/mol, which is three times higher than that of the reference RNR12 peptide (–2.29±0.25 kcal/mol) and formed a stable complex with imidacloprid. Furthermore, the dissociation constants (Kd), calculated from the binding free energy, showed that WQA32 (40 μM) has three orders of magnitude lower Kd than the reference RNR12 peptide (3.4 × 104 μM). Docking and RMSD scores showed that the WQA34 peptide is potentially selective to the target imidacloprid with respect to acetamiprid and clothianidin. Therefore, this peptide can be used in wet-lab experiments to prepare a biosensor to selectively detect imidacloprid.


Introduction
Imidacloprid (IMI) is also known as 1-(6-chloro-3-pyridylmethyl)-2-nitroimidazolidine, is a top-selling neonicotinoid insecticide, it functions in the target species as an agonist of nicotinic acetylcholine receptors (nAChRs) [1][2][3][4][5][6][7].It has been widely used for crop protection and animal health due to its excellent efficacy and selectivity to ravaging insects [8].Based on the following findings from in vitro comparative functional experiments, it has been hypothesized that imidacloprid may be less neurotoxic to mammals, including humans, than to the target insects [2][3][4][5][6]9].However, concerns over the safety of mammals have surfaced [1].According to Yamakuni et al., nAChRs may dysregulate endogenous ACh-elevated catecholamine synthesis in the chromaffin cells of the adrenal gland, ultimately resulting in an increase in adrenaline release from this endocrine organ [1].
Biosensors are currently sought for uses in different fields of applications such diagnostics, food safety and environmental monitoring because of their cost-effectiveness and amiability [10][11][12][13][14][15].Peptides have recently been re-examined and their uses as recognition element have been investigated [16][17][18][19][20][21].Peptides are known as stable, cost effective and easy to prepare through standard protocols and can be modified to fit specific applications [16,[21][22][23].In fact, peptides are suitable for detecting proteins [21,24], nucleic acids [25], bacteria [26] and chemicals compounds [18,21,27].In addition, peptides can be easily obtained by chemical synthesis methods and avoid the need for in-depth and laborious procedures such as in the case of antibodies [28].In particular, research activity in peptide-based biosensors is increasing.Based on the search for "peptides" and "biosensor" keywords in SciFinder engine (SciFinder, CAS), several thousands of matches have been found [28].
Detecting the presence of residual levels of pesticides such as imidacloprid, acetamiprid and clothianidin is necessary in many fields such as healthcare, environmental monitoring and food safety [29][30][31][32][33].Therefore, we examine the binding between IMI compound with several peptides, which were designed by choosing the interacted amino acids in the active site of the Lymnaea stagnalis Acetylcholine-Binding Protein Q55R mutant complex which is very similar structure to insect nicotinic acetylcholine receptors [34].These peptides will then be guiding models for creating peptide-based biosensors to detect the present of Imidacloprid compound.

Selection of the peptides
In our quest to identify concise peptides suitable for imidacloprid binding, we initiated our investigation using structural insights derived from the Lymnaea stagnalis Acetylcholine-Binding Protein Q55R mutant complex (PDB ID 3WTH) [34], accessible within the RCSB PDB database [35].This pentameric protein represents a pivotal target for neonicotinoid insecticides in the insect kingdom.The architecture of each monomer is intricately composed of a combination of beta sheets and helical structures interconnected by coil regions.The pesticides usually stacked Tyr185 with basic residues such as Lys35 in loop5 [34].
The starting protein is composed of five monomers in complex with five imidacloprid (S1 Fig) .Three peptides were selected from the binding areas of the different monomer/imidacloprid complexes using the "PLIP" webserver software [36], which enabled us to identify the possible binding sites between the protein and the ligand.The 3D forms of the peptides were retrieved using PEP-FOLD webserver software [37][38][39].

Peptides
Based on the results from the selection of peptides and the study of the different poses of interactions between the IMI and the three peptides, we found out that longer peptides have more stable complexes and higher affinity toward the ligand.Therefore, we intentionally provoked several modifications and mutations to create four new peptides having lengths from 21 to 34 residues.Also in this case, we used PEP-FOLD to generate the 3D structures of the peptides.More information about the modifications and mutations are discussed in the Results and Discussion section.

Molecular docking
To evaluate the affinity between these peptides and imidacloprid, we performed a molecular docking simulation using AutoDock VINA (v1.2.3) [40], with the peptides serving as the receptors.To ensure reproducibility, the entire peptide was positioned within a defined grid box, which had dimensions of 66 × 40 × 40 Å 3 and the hydrogens were added to the polar atoms of the peptides.Imidacloprid was utilized as the ligand, and it's important to note that all bonds in the ligand were specified as rotatable, allowing for flexible conformational sampling throughout the docking process.These specific parameters and conditions were employed to substantiate our methods and make them accessible for replication by others.The AutoDock Tools (v1.5.7) [41] was used to prepare the input files.The pH of the systems is 7. PyMOL [42] and BIOVIA Discovery Studio Visualizer [43] were used to analyze the established interactions between the peptides and the ligand.

MD simulations
To confirm the stability of the selected poses, MD simulations were carried out using NAMD software [44] for a duration of 480 ns.Beforehand, the designed peptide-IMI complexes were minimized and equilibrated using NAMD software for 10,000 cycles.Detailed MD simulations using the complex structures were conducted with the CHARMM36 forcefield [45].The complexes were solvated in a cubic box with keeping 10 Å between the complex and the box edge.NaCl ions were added to neutralize the system charge.For all simulations, temperature was set at 310 K and the pressure was set to 1 bar to closely mimic the general wet-lab experimental conditions.Subsequently, the fully temperature and pressure equilibrated systems were used as the initial configurations for the MD production.All simulations were conducted using a 2-fs time step.
Root-Mean-Square-Deviation (RMSD) was calculated with tools included in VMD [46] to check the stability of the studied systems.MM-PBSA method [47], implemented in the VMD CaFE plugin [48], was used to calculate the average of free energy from 3 different trajectories between 440 ns and 480 ns simulation time, after the convergence of the systems [47,49,50].
The binding free energy (DG slvd binding ) can be decomposed into the relative free energy of the solvated receptor-ligand complex (DG Cplx sol ) and the separated, solvated ligand (DG Lig sol ) and receptor (DG Rec sol ) as shown in Eq 1: This free energy (DG s ) can be further decomposed into four main contributions from the van der Waals (DE vdW ), the electrostatic interaction (DE elec ), polar solvation (DG pol ) and nonpolar solvation (DG nonpol ) as shown in Eq 2: The dissociation constants were calculated from the values of the binding free energy for the different peptide/IMI complexes using the online available applet (URL: https://protsim.github.io/protsim)[51].The dissociation constant (K d ) is linked to the free binding energy by the following equation: Where R is the gas constant, T is the working temperature K a and K d are respectively the association and the dissociation constants.However, based on the information provided in the article [51], we used the method of orthogonal distance regression to calculate the dissociation constant (K d ) and did not use the theoretical method mentioned in Eq 3. We followed the instructions in the article to input the raw data and leave the calculation method as the default value to obtain the values of K d , K a , and ΔG.We cited the article to acknowledge the authors and followed the terms of the license for any modifications made to the source code.

Selection of the primary peptides
The analysis of the different complexes between the monomers and imidacloprid showed that the existence of 7 amino acids (143A, 185A, 192A, 53B, 104B, 112B and 114B from the A and B chains) interacting with the target pesticide (S2 Fig) .Based on these findings, we selected the sequence from A185 to A192 to build the part I (named: YSP09) and the other from 104B to 114B to build the part II (named: DRM12).A third peptide (named WQW13), consisting of 13 amino acids located between TRP53 and TRP65 from chain C of the protein, was also selected since it interacts with imidacloprid, as confirmed by the PLIP analysis.These results suggest that the interaction between TRP53 and imidacloprid is facilitated by a π-stacking interaction between the imidacloprid molecule and the side chain of tryptophan residue at position 53 in the protein chain.Furthermore, the molecular docking results show that these 3 peptides have good binding affinity scores ranging between -3.17 and -3.80 kcal/mol (entries 1 to 3 from Table 1).

Improving the performances of the peptides
To improve the interaction between the peptide sequence and IMI, we combined three shorter peptides (entries 1 to 3 from Table 1) to obtain four longer peptides (entries 4 to 7 from Table 1).YSP09 was linked DRM12 from its N-terminal position to create a 21-amino acid sequence (YSM21), which was tested for interactions with IMI.For (PSM22) peptide, we started from invert order of WQW13 and then we modified it by adding six (PSHSSD) and three (LYM) amino acids in the C-terminal and N-terminal positions, respectively, to build a new peptide with 22-AAs sequence.PSW31 was built from YSM21, and the sequences were modified by adding the first five amino acids of PSM22 (PSHSS) in C-terminal and the first five amino acids of WQW13 in invert order (TTRQW) in the N-terminal position to form a sequence with 31 amino acids.The last peptide (WQA34) is obtained by adding the first six and the last seven amino acids from (WQW13) to C-terminal and N-terminal positions of a Peptide code is generated from the first two amino acids from the right followed by the last amino acid and total number of sequences.b This peptide is studied in literature [27] and is run for sake of comparison. https://doi.org/10.1371/journal.pone.0295619.t001 peptide YSM21, respectively.The 3D structures of the entries 4 to 7 with IMI are shown in Fig 1 .After that, the stability, and the affinity of the four combined peptides toward IMI was assessed using molecular docking.

Molecular docking
From From the summarized data in Table 1, several preliminary conclusions can be drawn: i) the docking scores of the short peptides are lower than those found for the longer peptides but remain in the same magnitude of RNR12 (-3.80 kcal/mol for WQW13 vs. -3.90kcal/mol for RNR12), selected by using modified phage display library [27,[52][53][54], RNR12 is a highly specific oligopeptide sequences for the recognition imidacloprid [27] and ii) the docking scores increased after the lengthening and the combination of the 3 primary selected peptides.Indeed, YSM21 has the lowest score among the four new peptides (-4.30kcal/mol), but it is higher than the best score found for the short WQW13 and reference RNR12 peptides.Thus, the increasing of sequence length of the peptides increases the binding energy, which is probably due to a decrease of the peptide flexibility, but these conclusions should be confirmed by more accurate MD simulations.The 3D structure of RNR12 and the interactions from molecular docking results are shown in S3 Fig, it shows that RNR12 forms two interactions with imidacloprid resulting from the interaction with ARG3 and LEU7.

MD simulations
Thermodynamic analysis.Further work has been done only on the longer peptides (YSM21, PSM22, PSW31 and WQA34) and we compared their results to that obtained with the reference RNR12 peptide.All the data, gathered in Table 2, indicate that the corresponding molecular systems are energetically favorable since they all showed negative free energy values.
More in-depth analysis shows that the reference RNR12 peptide performed well in comparison to the results from docking.Indeed, the RNR12 has an average free energy of -2.29±0.25 kcal/mol, which is higher than that obtained with longer PSW31 peptide and slightly lower than that obtained with PSM22.YSM21 and WQA34 resulted in average free energy values higher than that obtained with RNR12.YSM21 and WQA34 have average free energy values higher by 154.6 ±3.2% and 281.2 ±1.1% in respect to the value obtained for RNR12.Potentially these two peptides can serve for a better recognition of the target molecule than the reference peptide.Furthermore, energy decomposition showed that the main contribution to the free energy arises from van der Waals and the electrostatic interactions, and nonpolar free energy.The polar free energy contributed unfavorably to the interaction since all the values are positive but remain weak compared to the other contributions.For instance, the favorable contribution is 144.6% of the interaction energy while the unfavorable contribution is 40.5%.The small difference comes from two weak terms which are a favorable entropic one and unfavorable dispersion free energy term, neglected in Eq 2.
From the data gathered in Table 2, we can see that the YSM21, PSM22, PSW31 peptides have high dissociation constants (1.4 × 10 4 to 3.4 × 10 4 μM), which are in the same order of magnitude of the reference peptide.This denotes that these complexes are not very stable because their relatively low average free binding energies.Peptide WQA34 has a dissociation constant of 40 μM, which is 3 orders of magnitude lower than the reference RNR12 peptide, making this in-silico designed peptide a viable candidate for the recognition of imidacloprid and building sensors for it.
RMSD analysis.The RMSD analysis indicates that there are significant changes in the IMI-peptide interactions between the peptide-free and associated states.Specifically, the RMSD values of PSW31 and WQA34 show a substantial decrease of 50% and 74%, respectively, indicating a more stable interaction with the peptide in the associated state Fig 3A .The reference RNR12 peptide displayed a RMSD almost superposable to that of WQA34, showing that both peptides are stable.On the other hand, YSM21 exhibits only minor variations in RMSD that do not exceed 10%, suggesting that the peptide-IMI interaction is relatively stable.From Fig 3B, we can see that a large fluctuation of the RMSD occurred denoting the ligand keeps changing confirmation during the overall simulation time.PSM22 displayed a significant decrease in the average RMSD value of 30% in the associated state compared to the free state, indicating that the complex formation has a stabilizing effect on the peptide conformation.
Overall, these results provide valuable insights into the IMI-peptide interactions and their dynamics, which can be useful in understanding the underlying mechanisms of biological processes.The time-dependent RMSD analysis reveals that the peptides PSW31 and WQA34 reach stability after 100 ns of simulation, while YSM21 and PSM22 showed RMSD fluctuations till the end of the simulation time, although they are a large number of contacts (09) with the target, as shown in S4 Fig.This difference in stability can be attributed to the length of the peptide sequences, which increased from 21 to 34 in the case of PSW31 and WQA34.The longer peptide sequences may provide more stabilizing interactions with the pesticide molecule and reduce the fluctuation in the IMI-peptide complex structure, leading to a faster stabilization.Conversely, the shorter peptide sequences of YSM21 and PSM22 may not provide enough stabilizing interactions, resulting in higher structural fluctuations and slower stabilization.These findings highlight the importance of considering the length and sequence of peptides in studying their interactions with proteins and the dynamic behavior of imidacloprid-peptide complexes over time.

Selectivity
We also examined the selectivity of the lead peptide (i.e.WQA34) by comparing it with two other pesticides.The first one is acetamiprid (ACE) [30] and the second one is clothianidin (CLT) [7,9], which are two widely used neonicotinoid insecticides from the same family Fig 4A [9,30].As shown in Fig 4B , The MM-PBSA results of IMI, ACE and CLT were respectively -5.88±0.38 kcal/mol, -2.00±0.48kcal/mol and -1.02±0.55kcal/mol, demonstrating a higher affinity of WQA34 to the target pesticide although there is structural similarity between them.This is confirmed by the RMSD values (IMI: 0.3, ACE: 0.9, CLT: 1.9).Indeed, we witnessed a higher stabilization of the complex formed with IMI in respect to those formed with ACE or CLT.The RSMD score of CLT is approximately six times higher than that found for IMI and its docking score is lower by 66%.This suggests that the selected peptide is potentially selective to IMI in the following order IMI>ACE>CLT.

Conclusion
Starting from nicotinic receptor-imidacloprid complex available in the PDB databank, we selected three primary short peptides that can complex the target molecule with energies close to that obtained from a reference RNR12 peptide, selected in wet-lab experiment using Phage display.The combination of these peptides allowed their lengthening and increasing of the overall stability.The latter showed higher affinity to the target as showed by the MM-PBSA method.In particular, the WQA34 peptide has a threefold higher affinity and 850 times lower dissociation constant with the target compared to the reference peptide and showed that it formed a stable complex with imidacloprid as demonstrated by MD simulations.Furthermore, the peptide showed to be selective to the target pesticide against two other members from the Commencing with the nicotinic receptor-imidacloprid complex found in the PDB database, three initial short peptides were identified capable of forming complexes with the target molecule, exhibiting energies closely aligned with those of a reference RNR12 peptide, which was selected through wet-lab experiments using Phage display.Combining the peptides led to the lengthening of the peptide sequences and enhancement in overall stability.In particular, the WQA34 peptide exhibited a threefold increase in affinity and an 850-fold reduction in dissociation constant when compared to the reference peptide.Molecular dynamics simulations provided evidences that WQA34 formed a robust and stable complex with imidacloprid.Moreover, these investigations revealed the peptide's selectivity towards the target pesticide when compared to two other members of the same chemical family.Overall, these results indicate that WQA34 holds significant potential for application in wet-lab experiments, enabling the selective recognition of the target pesticide.However, it is essential to exercise caution when interpreting these findings, considering the inherent approximations associated with MD simulations and the potential influence of pH on the calculated values of binding free energy and dissociation constants.

Fig 1 .
Fig 1. 3D structures for IMI and the four peptides YSM21, PSM22, PSW31, and WQA34.In the middle is the structure of the IMI compound surrounded by four peptides YSM21, PSM22, PSW31, and WQA34 which have been obtained from the combination of the initial ones.https://doi.org/10.1371/journal.pone.0295619.g001

Fig 2 .
Fig 2. Interactions between IMI and a) YSM21, b) PSM22, c) PSW31 and d) WQA34 form the molecular docking results.The non-covalent interactions established between the IMI compound and the four peptides as the results of the molecular docking generated by BIOVIA Discovery Studio Visualizer software.https://doi.org/10.1371/journal.pone.0295619.g002

Fig 3 .
Fig 3.The molecular dynamic results: A) the root mean square deviation (RMSD) for the free and the associated peptides, B) The root mean square deviation (RMSD) for the peptides vs time.(The root mean square deviation (RMSD) in Å for the five peptides YSM21, PSM22, PSW31, WQA34 and RNR1: the free peptides (in red color) and the associated peptides (in orange color); B) The root mean square deviation (RMSD) for these five peptides vs time in nano second).https://doi.org/10.1371/journal.pone.0295619.g003 Supporting informationS1Fig.The structural of the Lymnaea stagnalis Acetylcholine-Binding Protein Q55R mutant complex (PDB ID 3WTH).(DOCX)

Fig 4 .
Fig 4. A) Chemical structure of the three neonicotinoid pesticides and B) Plots of root mean square deviation (RMSD) and the average free energy determined for the selectivity of WQA34 peptide (the error bars are obtained from three different trajectories).Figure (A) shows the 2D chemical structure of the three neonicotinoid pesticides IMI: imidacloprid, ACE: acetamiprid, CLT: Clothianidin and Figure (B) shows the plots of root mean square deviation (RMSD) in Å (grey color) and MM-PBSA results in kcal/mol with the error bars (green color) determined for the selectivity of WQA34 peptide.
Fig 4. A) Chemical structure of the three neonicotinoid pesticides and B) Plots of root mean square deviation (RMSD) and the average free energy determined for the selectivity of WQA34 peptide (the error bars are obtained from three different trajectories).Figure (A) shows the 2D chemical structure of the three neonicotinoid pesticides IMI: imidacloprid, ACE: acetamiprid, CLT: Clothianidin and Figure (B) shows the plots of root mean square deviation (RMSD) in Å (grey color) and MM-PBSA results in kcal/mol with the error bars (green color) determined for the selectivity of WQA34 peptide.https://doi.org/10.1371/journal.pone.0295619.g004

Table 1
On the other hand, PSM22 has more affinity to the target since it docking score is 35% higher (-4.61 kcal/mol) resulting from several non-covalent interactions with imidacloprid, within 4 Å.We can notice two conventional HBs, one from GLN18 (2.61 Å) and one from ALA8 (4.00 Å) and two carbon HBs from TRP15 (2.40 Å) and from GLN18 (2.38 Å).The binding energy of peptide PSW31 is even higher than that found with PSM22 (-4.80 kcal /mol).The receptor established two carbon short hydrogen interactions (<2.70 Å) with SER5 (2.66 Å) and with SER7 (2.60 Å).WQA34 has the highest score docking among the four peptides under the study (-5.30kcal/mol), although it shows only one HB with the target from TRP6 (3.91 Å) Fig 2.

Table 2 . The average free energies result in kcal/mol for the solvated systems, details of the energy contributions for the different complexes and the dissociation constant a .
The average free energy of solvated receptor-ligand binding (DG binding;slvd ) and the four main contributions from the van der Waals (DE vdW ), the electrostatic interaction (DE elec ), polar solvation (DG pol ) and nonpolar solvation (DG nonpol , equilibrium dissociation constant (K d ).